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ABSTRACT 


It is the purpose of this paper to develop a useful mathematical 
model of ASW aircraft availability. The increasing emphasis of systems 
studies dictates the use of accurate and representative models of the ASW 
systems. At present, many studies are using essentially the same models 
developed during World War II. This paper is an attempt to make use of 
advanced theory in a more powerful and flexible model and to make the use 
of the model practical and verifiable. 

The writer adapted the time homogeneous bivariate model as de- 
veloped by F. C. Collins. This is a discrete time Markov process with 
a stochastic matrix of transition probabilities wherein the maintenance 
process is modeled as a pulsed input multiple server queue, 

The model was programmed in FORTRAN 63 on the CDC 1604 and 
then modified to allow for variability.in the input parameters. Other 
modifications include an increase in the size of the model to accommodate 
a l6-aircraft squadron, the largest ASW squadron at present, and an 


explicit form solution to the maintenance queueing equations. 
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Symbol 


а/с 


X, (t) 
X, (t) 


A(t) 


N(t) 


TABLE OF SYMBOLS 


Definition or Meaning 


aircraft 
mean repair rate of aircraft 
mean accident rate 


the set of all possible outcomes; the probability 
description space 


possible outcome(s) or event(s) 
a conditional probability that at time n + 1 the 
outcome or state is j given that at time n the 


state is 1 


the number of a/c flying at time t, which did not 
fly the previous cycle 


the number of a/c in the maintenance queue at 
time t | 


the number of a/c desired on station at time t 


the total number of a/c of type considered at 
time t 


the time interval from the launch to recovery at 
the start of the cycle 


the probability distribution over all possible 
states at initial time to 


the matrix of transition probabilities at time t 





Symbol 


P (o, i) (8, j) 


Ygh 


Definition or Meaning 


ΤΠ or the E matrix; the probability that 
х, aa el X, 


that A; = @ апа A, ο ο [be start oí the cycle 


= j atthe end of a cycle, given 


the probability given f ready a/c, g are launched, 
and h enter maintenance 


paopabllity of entering maintenance just before, 
during, or immediately after launch 


the probability of equipment failure during flight 
requiring maintenance when recovered by the 
earrier 


the probability that of a a/c flying m will enter 
maintenance upon recovery 


the number of independent identical maintenance 
repair stations or ''spots'' 


the probability that i - j a/c are repaired in 
time interval t 





l. INTRODUCTION 


The threat to freedom of the seas posed by the vast Soviet sub- 
marine fleet is perhaps the most thorny problem facing the U. S. Navy 
today. Two world wars have produced Pyrrhic victories over limited 
submarine fleets. During the Second World War operations analysis 
was born into the Navy to aid in the defeat of the German submarine. 
The classic antisubmarine warfare (ASW) analyses and models developed 
by Morse [2] and Koopmans [3] are still being used today, over two 
decades later, in most of the ASW study efforts for the Navy. 

These early ASW analyses assumed a given level of search effort 
available and directly evaluated the probability that an ASW subsystem 
could detect and/or kill a submarine. This assumption is not only 
logical to make the problem tractable, but also practical since no 
immediate changes in ASW force levels could be expected. Moreover, 
the studies were conducted during the war, not before it started. It is 
the purpose of this paper to present a probabilistic model to describe 
the available effort. Such a model can be used to sharpen the estimates 
of the effectiveness of an ASW subsystem and to study the charac- 
teristics of the associated support system. 

Naturally, the current study plays an important but limited role in 
the overall problem of designing an entire ASW system. The diffi- 


culties involved in such a specification are legion. First and foremost 
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is the quantification of the ASW mission in denying the enemy the 
effective use of his submarines. Currently, the probability of detecting 
and/or killing submarines is used as the measure of effectiveness of 
the mission, and it appears that a more encompassing one has not been 
developed. Second, the specification of an ASW force level to counter 

a given threat has many inherent subjective elements. These are due to 
the existing historical bias in predicting the conduct of a future ASW war 
with an enemy, particularly one who has never before used a large 
submarine force in its military operations. The reader can imagine 
why merely defining terms such as 'threat!'and ''effective counter"! 
becomes quite difficult. 

Thus there is a needto investigate the levels of search efort 
specified. This may require acceptable models to measure the availa- 
bility of effort, its effectiveness, and determine the logistic support 
required for any level of available effort. Specifically, the ASW sub- 
system to be modeled is the carrier-based aircraft, although the model 
is adaptable to other systems. 

The method of investigating the demand for ASW carrier a/c will 
assume that the desired number of a/c on station is known as an input 
parameter. The support required to achieve this measure of available 
effort depends upon maintenance space, manpower, and supply. 
Generally, we shall consider how an ASW carrier supports this number 


of a/c on station with the present or proposed number of a/c embarked 
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onthecarrier. The. parametric input can be subjected to sensitivity 
analyses. 

The operational commander of the ASW force launches the desired 
number of a/c on station to screen, search, or actively prosecute a 
submarine contact. Each a/c is relieved on station. Each such relief 
requires the launching of another a/c prior to the recovery of the initial 
a/c. The returning a/c must receive varying degrees of maintenance 
and requires refueling and rearming. This cycle continues until the 
mission is completed. Loss of a/c due to accidents, insufficient supply, 
and lack of repair capability cause deviations in this procedure. Naval 
operations involve the interaction of many quantities which are random 
іп nature. Not all can be considered in a tractable mathematical model. 
Some quantities which are important are omitted. One example is the 
length of each cycle time, which is assumed to be a constant value. 
Including variables of this nature incurs unnecessary mathematical 
complication. It is hoped that adequacy of the model can be measured 
by using fleet data available from the Fleet ASW Data Analysis Program 
(КАРАР). 

Collins [5] describes a bivariate Markov model for airborne early 
warning (AEW) and combat air patrol (CAP) jet a/c operating in an 
attack carrier force. This model is used to evaluate the probability of 
maintaining a fixed requirement of a/c on station as a measure of ef- 


fectiveness of the system. It has subsequently been used in a larger 
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attack force study for the Navy. The model computes the probabilities 

of the number of a/c on station and in or awaiting maintenance at any 

given launch period. The comparable ASW problem differs in the following 
aspects: 

1. Type, range, and speed of a/c; 

2. The variable number of a/c required for mission; 

ic Attrition due to accidents and supply failures; 

4. The greater number of ASW a/c. 

It was decided to use the Collins' model with appropriate modification. 
For immediate reference, the mathematical content of the model will 
be repeated herein. 

In order to incorporate these modifications, ít was necessary to 
spend some time reprogramming on the CDC 1604 digital computer in 
FORTRAN 63, the CDC version of the IBM FORTRAN IV. The original 
program was not readily avallable and was written in an early assembler 
language. Moreover, the numerical analysis was not sufficiently sharp 
to handle the larger imput values. Also, double precision (two computer 
words instead of one) arithmetic was required in one subroutine for an 
accurate explicit solution to the maintenance queueing equations (see 
Appendix 1). This effected a 50% decrease in the computer time 
required for developing a matrix of transition probabilities. 

Following this introduction, section 2 contains a brief description of 


the operational problems involved and the assumptions made. A brief 


description of Markov chains and the mathematical model are presented 
in section 3. The details for computing the matrix of transition proba- 
bilities are given in section 4. General employment of the model 
follows. The appendices include the solution mentioned on the pre- 


ceding page, a logical flow diagram ofthe program, a copy of the 


program, and some sample results. 
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2. ASSUMPTIONS 


The real-world employment of carrier a/c is cyclic in nature, and 
the present state of any given a/c (i.e., flying, in or awaiting mainte- 
nance) depends largely on what the previous state was. This fact 
suggests that a Markovian assumption can logically be made for the 
a/c transition probabilities. In the search phase, a/c may or may not 
relieve on station; but, in any part of the contact investigation phase, 
relief on station will be made. To insure full screening and mission 
coverage, a/c will relieve on station. 

The question of resupply during an operation depends primarily on 
the availability of carrier on-board delivery (COD). This depends on 
the geographical location and the mission (convoy protection, strike- 
force protection, hunter-killer operation, etc.). In practice, resupply 
is not anticipated within a week's period, and around-the-clock oper- 
ations have continued for two weeks without resupply. 

Standard maintenance procedures aboard carriers preclude major 
maintenance on the flight deck. It will be assumed that sufficient notice 
is given so that all major 120-hour checks will be completed prior to 
the operation. This assumption can be modified with an appropriate 
adjustment in the mean repair rate. The concept of maintenance crews 
assigned to hangar deck areas (''spots''), as developed by Collins [3], 
will be used. Each crew will be capable of all types of maintenance 


2 OL 
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and will operate independently at the identical mean repair rate \. 
The number of spots is determined by the average number of such 
crews available to work continuously around the clock on a watch basis. 

The state of each a/c is assumed to be statistically 1 of 
that of others, and the launching and landing transition probabilities 
will be developed on the basis of independent Bernoulli trials. The 
parameters can be determined using the maximum likelihood estimators. 
The range of the number of a/c desired on station at any given cycle will 
be the user. The number to be launched at any time is assumed 
equally likely within this range. This input parameter is a function of 
the estimated ETE density (i.e., expected contact rate). The 
lower limit will be set at the number of a/c desired on station in the 
search (screening) phase, and the upper limit is set at the maximum 
practicable number of a/c to be launched during a multiple-contact 
phase. 

Briefly, the assumptions are: 

l. a/c will be relieved on station. 

2. Any desired length of operation can be set as an input. 

3. Major 120-hour checks will be completed prior to the operation. 

4. No resupply to the carrier is available. 

5. The launch-to-launch cycle for all ASW a/c is four hours. 

6. Minor maintenance, refueling, and rearming only can be 


performed on the flight deck. 


ju 


m" 


Each maintenance spot is characterized with an independent 
exponential repair time with mean repair rate of A for 
around-the-clock operations. 

The number of a/c lost due to attrition is a Poisson random 
variable for each cycle period with parameter ha (8/ς 
accident/flying hours for a/c type). 

Any a/c lost by accident will not be returned to service due to 
either (a) physical loss at sea, or (b) insufficient mainte- 
nance capability aboard ship and lack of major parts, 

The number of a/c launched for each cycle is uniformly 
distributed between the upper and lower limits determined by 


the user. 
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3. MODEL DESCRIPTION 


sene Theory 

A stochastic or random process is a collection of random variables 
indexed on some set T, (X(t), te T). In this case, time is the indexing 
set, and the Markovian assumption states that the future state of the 
process depends only on the state at the present time and not on its 
past history. Due to the cyclic nature of our problem, it is possible 
τυ ΕΠ IND) ) Using the cycle time from launch 
to launch as the steps of unit time in a discrete Markov chain. It is 
assumed that the reader is familiar with the notion of a random variable 
as a function defined on a sample description space (S) on which the 
family of events or outcomes (E) ips probability function can be 
defined [4]. 

A discrete time Markov chain is described by a sequence of dis- 
crete valued random variables and is determined when the one-step 
transition probabilities of the state variables are specified, i.e., a 
conditional transition probability of a transition at time n for each pair 
of i j=0,1,...,m (me being the number of states in the process) 
must be given. 

P; Ut nm t. 1) = P [X (m + 1) = j| X(n) =i] 
If the transition probability functions depend only on the time difference, 


we have time homogeneity 
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Е, 1) = peo) = eer. 
ΠΠ ο ο ο. 


The initial state of the system must be given either as a specific state 
or randomly as a probability distribution function over the possible 
states, 

Ше Pi; (transition probabilities) are arranged in matrix form and 
satisfy: 

| во form 0. DE 


m 
2. 5 p.. = 1, ie., the rows of the transition matrix sum to 1 


jo ^? 


for all i for the states within the description space [4]. 
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3.2 The Model 
In order to establish the finite set of states (E) for the model, we 


shall consider two random variables defined as follows: 


x (t) = The number of a/c flying at time t not having flown in 
the previous launch-to-launch interval. 
X.) (t) = The number of a/c in or awaiting maintenance at timet. 


т - 
$ + 


Now, we will consider the vector X(t) = LX, (t), X, (t) ] as a pair of 


random variables and thereby have a bivariate stochastic process with 


M: 


the possible states ranging from (0, 0) to (A, N). 
ος Xi (1) < А = No. of a/c desired on station, and 


ος X, (t) < N No. of a/c of given type aboard carrier., 


We will define an operating cycle as an interval unit of time. Process 
observations of X(t) will be made at successive unit interval launch 
times. To develop the Pi; elements, consider a given time t for launching 
until A aircraft are flying or until the supply of ready a/c is depleted. 
Those a/c failing the launch enter the maintenance state at this idealized 
point in time t (the total launching time required is much less than the 
total cycle time. At some time T, less than the launch-to-launch unit 
time interval, the a/c which were relieved on station return and land at 
the idealized point in time t+ T. Some of these a/c will require mainte- 
nance and enter the maintenance queue. Those requiring only refueling 
and preflight inspection will enter a ready status to be tested for the next 


launch. 
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During the unit time interval, maintenance will be performed on 
those a/c in the not-ready status, and a certain number of aircraft will 
be repaired according to assumption 7. 

In summary, we start the system in some initial state (such as 
(0, 0) with no a/c flying or in maintenance) or start with a probability 
distribution Q(t) over the states, E, at time = . We launch, recover, 
and repair a/c in the unit interval and repeat the process over each 
succeeding unit time interval until the end of the operating period. 
Knowing the transition probabilities within the unit time interval, we can 
develop the elements of the transition matrix, P, or M і), (В, 3} : 
These are the probabilities of going from the state of q a/c flying and 
i a/c in maintenance to § a/c flying and j a/c in maintenance over the 
unit time interval. 

It was assumed in section 2 that A, the number of a/c to be launched, 
and N, the total number of a/c on board, are random variables, whereas 
they have been treated as constants so far in the development, To be 
analytically correct in including this feature, one πο develop the 
appropriate quadrivariate process. Such a development leads to too 
large a state space and the author chose to include these effects by using 
a Monte Carlo simulation technique. That is, at the beginning of each 
cycle, a random mechanism is used to determine the values on À and N. 

The probability of losing an a/c or changing the desired number to 


be launched is determined from the specified distributions at the beginning 


dl 


of each unit interval, and the resulting P matrix containing the 
БТ i), (B, j) is then recomputed. The probability distribution Q(t) 
over the states at any time t may be determined by the appropriate 
number of successive iterations of the Q vector times the P matrix, i.e., 
Q(t) = P[X, (ty) = 6. X, (t) το 

The probability of maintaining q a/c on station over any given period 
of operation may be obtained at any unit time t (i.e., the beginning of 
the next cycle) by summing out the appropriate maintenance state proba- 


bilities. Thus, P(o a/c are flying at time t ) = 


EL (X, (t) 2 o) = È Pr (X, (t) = 0, X,(t) =i) 
t=O 


A mathematical comment appears to be in order. In the case of 
fixed A and N, the states of the Markov chain are positive recurrent; and 
steady-state probabilities can be лын for the entire state space. In the 
case of decreasing N due to a/c attrition, this is not true; and (0, 0) 
becomes an absorbing state as time (t) goes to infinity. This latter 
consideration is not a realistic one for the operational period envisioned. 
Therefore, it is mathematically more feasible to use the former chain in 


conjunction with the Monte Carlo technique. 
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4. 


DEVELOPMENT OF THE TRANSITION MATRIX 


Perhaps the simplest way to view this development is to note the 


various transition probabilities incorporated in one-unit time cycle 


defined as follows: 


(1) 


(2) 


Veen = the launching transition probabilities at time t. This is the 

probability of taking f ready a/c, launching g successfully, and sending 

h into maintenance. Each a/c to be launched is considered a 

Bernoulli trial with probability of failure of p 7 which is estimable 

and subject to sensitivity analysis. The values of M are: 

a. O0 ifg» A, since only A a/c are desired; 

b. 0 ifg+h> f; it is impossible to launch and send into mainte- 
nance more a/c than are available; 

c. 0 ifg<A, g+h<f; launching continues until A a/c are flying 
or until all f are used up; 

d. К (1 - Р)” D y 8 if g < À, g+ h = f, standard binomial 

when all a/c in the ready state are used up but the A a/c are 

not launched; 

Er Ἡ (1 - p) ° {с if g=A, g+th>f, standard negative 


binomial for g successes in g +h - 1 trials. 


Hn (m) = the landing transition probabilities which occur at time 


tr T. We must consider the probability that if there are a/c flying 


at time t then m a/c will enter maintenance at recovery time t+ T. 
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П, (m) will equal a standard binomial where p = the probability of 


eguipment failure in flight: 


| α o -m, m 
= 1 - ӘЛЕ 2220. 
1 (m) = (a) (l -P (p , m=0,1 α 
(3) Fx (τ) - {πε maintenance transition probabilities, i.e., the proba- 
bility of repairing (i - j) a/c in time t. Two maintenance periods 


occur: the first starting at time t and ending at time t + T, the 
second starting at time t + T and ending at the end of the cycle, 

(t+ 1). Under assumption 7, the pulsed input, multiple exponential 
server queue is developed with D maintenance ''spots'' or servers 
each with identical, independent service rates, A. For each server, 
then, the probability of remaining occupied (given the server is 


busy) in time T=e ыт . The probability of becoming free (i.e., 


repairing апа/с)-1-е hah The resulting queueing equations are: 


Аза αι 


” а (9 / dt 


-лАР, Q(t) * (n+1) AP, 


< . 
i με U) for 0 < n < D; 


B. dE 


; „(t)/ dt 


- DAP. (t + DAP. (t) for nm > D . 
1, n 1, 1 


mt 

Three ranges of 1 (initial queue state), j (final queue state), and 

D become significant: 

a. When j<i< D, then not all spots are busy since there are 
fewer a/c in maintenance than spots. Each spot works inde- 


pendently; therefore, the solution to A is the binomial: 


ag 673 


p (0 = 1} =e) pee 
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When D < j <i, then all spots are occupied throughout the total 


service time, and the closed form solution to B is the Poisson: 


(my πας 
|. (DA t) e 
de 0-0: 


When j< D< i, then all spots are busy at the beginning of the 
service period, and some spots become idle during the service 
period. The explicit form solution of equation A is found using 


moment generating function transformation: 


ln (9) (D) (G= -) (ο η 


n=) 
— — 1-D-k 
чт } 


(The derivation of this solution is discussed in Appendix I.) 





The figure on the following page will show the relationships of these 


transition probabilities within the unit time interval. 
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TRANSITION PROBABILITIES WITHIN THE UNIT CYCLE 


FIGURE 1 
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In order to develop each transition probability over the total unit 
time interval, we must consider all events taking place within the 
interval. Thus, to obtain the probability of going from q a/c flying and 
i a/c in maintenance to 8 a/c flying and j a/c in maintenance, we start 
at the state (w, i) at time t. At this time, a/c are launched and some 
l a/c failing the launch enter maintenance. These i+ l in maintenance 
are then serviced until time t * T when some k a/c are still in the 
maintenance state. At time t+ T, ofthe œ a/c previously flying, some 
m enter maintenance and (o - m) enter the ready pool. Maintenance is 
continued on the (k * m) a/c for the remainder of the cycle (1 - T), until 
the end of the unit time interval when j a/c remain in the maintenance 


state, In functional form: 


N-o-i itl o 
Σ Σ 2 


Y И 
l=zo К=о m=0 N-q@-1, 8, 1 


Pia, i), (8, j) 


T) . (m) (1 SEES 


ΙΤ ' Pk+m, j 
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5. SUMMARY 


Representative values for the mean repair rate and the landing and 
launching failure rates produced results in agreement with the sensitivity 
analysis by Collins on these parameters in [5]. For failure proba- 
bilities less than . 5, and mean repair rate less than 12 hours, the effect 
of reducing the available maintenance time to 80% of the cycle time was 
negligible. Optimal loading and cycling policies can be determined for 
known values of these rates, 

The model affords the following checks: (1) the rows of each P 
matrix are summed as they are computed by the program; and (2) the 
probability distribution vector (QJ) is summed over the states, Each 
summation was within 10 Р of one a the computer model. 

The user may substitute any available distribution over the interval 
of a/c desired on station, In order to keep A fixed, enter the desired 
value as both upper and lower limit ( A - ALOLIM - LUPLIM). For 
fixed N, use a very small value for ALAM (such as 10 5 ). Subroutine 
KRAN is a uniform generator, using the half open interval (lower limit 1, 
upper limit + 2) and a starting pnm as inputs. KRAN outputs an integer 
in this interval. Subroutine DRAW was used to provide some intuitive 
grasp of the results. DRAW was used in binary card form and is not 


essential to the main program. (The indicated associated statements 


must be removed, however.) 
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The results of reasonable arbitrary parameter values, based on the 
author's experience, have shown that most of the probábilities concen- 
trate over a few states. Moreover, cómputation time increases rapidly 
as a function of N (no. of a/c), see Figure 2, This would indicate that a 
simple approximation to the model could be developed. One method 
presently being investigated to reduce computation time is to shrink the 
probability state space to include only those significant states and, thus, 
reduce the size of the transition matrix. Alternatively, the eigenvector, 
eigenvalue representation of the P matrix, might be used. 

Originally, it was hoped to utilize the data m the Fleet ASW Data 
Analysis Program (FADAP) to attempt a verification of the model with 
its real-world counterpart. The only method available at present for 
obtaining the necessary data is by direct observation or a program of 
data collection, as suggested by Collins [ 5]. 

Many fruitful areas of investigation exist: 

(1) Attrition has been simply modeled by the Poisson method. The two 
components of attrition, accidents and supply shortage, can be more 
accurately modeled and used to develop logistic schedules for 
maintenance and supply. One simple technique is to assume each 
component is Ж шде: a Poisson, and estimate a supply failure 
rate for AOCP η.» from past data. With these assumptions, 
the total attrition is Poisson, with the parameter equal to the sum 


oí the accident and supply failure rates, 
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(2) 


(3) 


(4) 


(5) 
(6) 


The model can be modified to make the number of maintenance spots 
available for any cycle a variable function of time, D(t). 

An investigation of the Markovian assumption validity as the cycle 
times become smaller and smaller. 

Development of a continuous time model, 

Modification of the model to simulate resupply by COD. 

A study of the distribution of submarine contacts to determine the 


validity of the uniform a/c demand assumption. 
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Time (min) 


PROGRAM ASSEMBLY AND COMPUTATION TIME 
FOR ONE TRANSITION MATRIX (P) AS A FUNCTION 
OF THE TOTAL NUMBER OF AIRCRAFT (N) 


FIGURE 2 
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Total Number of Aircraft (N) 
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APPENDIX I 


EXPLICIT SOLUTIONS OF THE 


MAINTENANCE QUEUEING EQUATIONS 


The queueing equations for the pulsed input queue are essentially 
the pure death process given in[1] and [4] as problems and developed 


by Collins in [5]. The equations are: 


dP. O 
e c 2. D 
A. T же) + τπτ for O<n< 
dP. , V) 
B. a to Τη + τ for n2 D 


where ED - SE and Е (Е) = 0 fori<j, since no input (arrivals) 
occur during the service time. 
Equation B is solved directly in closed form: 


(DA t) (n - 1) ve X Dt 


uu Να (n - i)! 


Now transforming the first equation (A) using the moment gener- 
ating function (MGF), 
р-1 , 
G(s,t)= E s P (t), 


n=0 


as outlined in [4] (Chapter 7), and its partial derivatives: 
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dG D-1 п / 





(1) PIN Ξ > s E (t) 
Ш-Н 
D-1 
(2) ΞΕ = Z ap (t) 
ds -0 n 


Where P (t) denotes the conditional probability Р. та) , by substituting 
п , 
(A) into (1), properly identifying the first summation with (2), and 


changing the second summation index to r - nt 1, we get: 





D 
- 1 
<= οι τοι πως Pa (t) e 
r= 0 
dc" dG E 
(3) $E - -A(5-1) X + ADs ΕΠΟΝ, 
since 
D 
r - 1 dG D - 1 
= rs P (t) ΠΕΙ + Ds PUO 


Next, replace the partial differential equation (3) with a system of 


ordinary differential equations using the Lagrangian auxiliary equations: 


The solution to the first equation (using the first two differentials) 
15; 


toss Ds e 
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and hence 


Επ. eS odi 
Or 
с, = dE ο.) 
The second equation is: (using first and third differentials) 
At D-1 
dz = -AàAD(C.e + 1) ΠΠ 
1 D 
Using the solution to (B) where m = 1 - Dto replace P (t) and 
| | : At D-1 
integrating, term wise, the binomial expansion of (€i e zd 
In spp | Р 
215 АШПШ) Σ T" Ge | зак 
m! 1 1 


j=0 
where the integral is evaluated as: 


k уре у 


> ї е aa + С 
m-k>v+wl k! 2 
k=0 (A(D - j)) 
Thus, 
D-1 m k m - k 
M t 
j=0 k=0 | = 
and the general solution is ¢ (С, : с.) , where ¢ is an arbitrary function 
and 
C - u(s, t, z) 
and 
С, vA SL NE 
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To get our particular solution, use the boundary conditions for G(s, t): 
Oor s=: 1, 


D-1 
GPO = 27 Р (0) 
п= 0 


= Pr [no. in maintenance at t is< D| iat t= 0] 





i- - À Dt n 
À Dt 
Gd t) 1 — Z = Mm Mo 
n! 1 
п=0 
where 
А = C =m 
m k 
E NN - С, = 2+е ^ D* Һр 
k=0 
so 
с. = z + y, (t) 
(2) fort 0, 
D-1 
G(s. 0) = 2 s P (0) = 0, since i> n> D 
n=0 
where 
u (s, 0, z) = (s - 1) 
D-1 : 
v (s. 0772) = C. =E L. (s - 1) (s) 
2 J D - j 
3=0 
Thus, 
D-1 i m 
S D j D ч 
GS о (7) с η е7 
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Substituting the general value for C, above: 


D-1 И : m 
D - At D 
G(s, t) = фы, у) = Z (jp) (s-1? e ^ (στ) 
| J DES 
Hs 
Du D cs (a Dt) * - ADt D p 
= > (т) (5 = 1) > k! e € > 
j=0 k=0 4 
Rearranging terms, i 
D-l | D-1 ү 
εν, ΙΙΙ С” (ы, Су е жа 
п=0 j=n 
& 
"D т Оро“ = se = | 
k=0 Don 
where P (t) = the coefficient of s” 
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KPEANMDUUA IT 


THE LOGICAL ГТО DIAGRAM OF THE COMTUTER PRCGRAM 





V A E 
eam / ES ») ρω e - — 

ο ο, ο ο οσα ICYCLE. JOYGLE 

y INPUTS | Вл DENEN LLOLTM, ETC.) 

CUN V 


9 


mu c't. А 


805 > DRA ТА) 
ç 


πα 


is Зав 


BE ves x TEIS EO 
dm NO CYCLE? 





(lower limit s A<upper limit) 












no қ там 


TUM -ᾱ-- Teen an 


—— m TA=TA-2 


y 


LLL yes TARTA. 
ее и ος. t=- :- ΤΑΞΤΑ-2 
| ё 


Ф ч TIA — w 9 


(1)----------------η SEL CANE y.) A — . —P 


* 
. 
— 





-T € E 
| COMPUTE "Т (1)=PRFI.(IAA, EL) | 


ee nn A 


ect ыл 







COMPUTE JIT) 





p13(t) 9 PTRE(I,T 


E or 
| | 
COMPUTE e, -GAMZA(IF2,IO,IH): | 


n 














— 





| | / 
ΘΕ... 


~ SET THE TRANGTITOV MATIX =0 
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801 


pr 


* 
p m — a Q = м ag cu dpa - s q quas сы 


| TENIS = Y 
| COMPUT. Hy ,x)( η) Ὁ) 1 
= 





< 


€ — m glass 


PIQUE EUM A Sm 


T 


= 





— 


= АА επ ποτ ο. 


A A A A A —À 


PRINT TIE p¿¿'s OF THE (IALPY,IX)TH ROW AND ROW SUM | 
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SUM THE BLIM! 


iul: eyes 
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m 
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SET THE QJ V;CTOR =0 
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А/ 


x 15 DS E | 


APPENDIX III 


THE COMPUTER PROGRAM 
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-COOP 5 sLANMAN 90/49/S/15/2S/F/45=54 5175 930000950 000 






-ВІМАВҮ, 56, 000 
(RELOCOM, 000 
-ЕТМ»Е. 000 

PROGRAM MARKOV 000 

C THIS PROGRAM IS A NONSTATIONARY BIVARIATE MARKOV CHAIN MODEL OF ASW A/C 000) 

C OPERATIONS. THE RANDOM VARIABLES ARE THE NUMBER OF A/C FLYING AT THE 000 

€. BEGINNING OF ANY GIVEN LAUNCH CYCLE. THE МАХІМОМ №. ОЕ А/С ALLOWED IN 0004 

C THE MODEL IS 16(МА). ТНЕ RANGE OF A/C TO BE LAUNCHED AT ANY GIVEN 0001 

SEMEINTERVAL 15 O TO 6 A/C. THE FOLLOWING INPUTS ARE REQUIRED, 0001 

C ID- THE NO. OF INDEPENDENT MAINTENANCE SPOTS 001 

С NA= TOTAL NO. OF A/C TYPE ON BOARD 001 

d TIMEZTIME FROM LAUNCH TO RFCOVERY/LAUNCH TO LAUNCH CYCLE TIME(HRS)001; 

(с FLAM=MEAN REPAIR TIME PER SPOT/LAUNCH TO LAUNCH CYCLE TIME (ΗΒ5)001. 

С PGAM= PROBABILITY OF A/C FATLING LAUNCH(MeLeESTe FROM PAST DATA) 001. 

@ P= PROBABILITY OF A/C FAILURE DURING FLIGHT REQUIRING MAINTENANCE 001! 

C AT LANDING (MeLe ESTIMATOR FROM PAST DATA) 001 

С GI- THE PROBABILITY DISTRIBUTION VECTOR OVER ALL POSSIBLE STATES 001: 

С (7 X 17 = 119) SUCH THAT THE SUM OF ALL QI(I) = ie THIS 0011 

С IS ESTIMATED BY THE USER AND INPUTTED BY USING A DATA STATEMENTOO1' 

С ICYCLE = NO. CYCLES DESIRED FOR OPERATION 00 2 

G JCYCLE = LAUNCH TO LAUNCH TIME(HRS)(TOT, TIME=ICYCLE x JCYCLEË) 002: 

C ALAM = ACCIDENT RATE FOR TYPE A/C (ACCIDENT/HOURS)y' 002: 

G ALOLIM = DESIRED LOWER LIMIT ON A 002; 

С AUPLIM = DESTRED UPPER LIMIT ON A 002: 
COMMON FLAM,TIME 002%} 
TYPE DOUBLE FLAM 0024 
COMMON PTRM, GAMMA PR ,PRFMA , ID 002: 
ЕЕ ЕТА ЕВС т : 0026 
DIMENSION PTRMCITS 1752) :GAMMA C179 7517) SPRFMA (7 97) 0027 
ОТМЕМ5ТОМ РЕ(119,7,17),01(119),0)(119) 002€ 
DIMENSION FIPLOT CITOJ ITI) 0025 

С ENTER DATARCARDS HERE | 003€ 
DATA((QI(I)>I=1>119)=22>16(.05)>102(,0)) | | 
NA=16 
А(АМ-.01 
ID = 8 

FLAMz3,0 
РСАМ=Р=.4 

ТҮҮ = 13421775 
ГІМЕ = .125 
ICYCLE=20 
MCYCEESA .... 
ALOLIM=4, 
AUPLIM=6, 

с END OF DATA CARDS 0031 
ALZALOLIM-1. $ AU=AUPLIM +25 | 0032 
UNITT=le 0033 
N=NA+1 0034 
ІАМАХ-7 0035 
IALAST=O 0036 
D=FLOATF(ID) 0037 
NLAST=NEWN=N ' | 66 38) 
КТ= 1 | 00391 

809 ITA=KRAN(AL»AU,IYY) 0040 


I (KIsEQ;1) 1139115 0041! 
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483 





50 


303 
300 
АТ 


FROM THIS NEXT STATEMENT TO NO. 


T1=-LOGF(.000000001 * RANF(-1))*2430258/ALAM 
IF (TI-TFLC) 13051315132 

T22-LOGF(4,000000001 * RANF(-1))*2430258/ALAM 
ШЕСТІҒТ2-ТҒІС) 230,;2315121 
T32-LOGF(4,000000001 + RANF(-1))*2.30258/ALAM 
IF (T1+T2+T3-TFLC) 331593319231 

NEWN=NLAST-3 $ GO T0113 

NEWN=NLAST-2 39 .σο T0113 

NEWN=NLAST $ GO 70113 

NEWN=NLAST- 1 

PRINT 8882, IA»NEWN 


IF(NEWN-IA) 15513913. 
TA=NEWN 

IF(IALAST) 11912911 
CONTINUE 


483 IS CONCERNED ONLY WITH THE GRAPH 
DO 482 1=1,12 

JT(1)=8H 

JT(1)=8HE(A/C) = 

JT(3)=8HSPOTS = 

JT(5)=8H T = 

JT(7)28HJ VS QJ 

JT(8)=8HVECTOR 

JT(9)=8H N = 

JT(11)=8H А = i 
DO 483 1=1,119 ς 
FI=]1 

FJPLOT(I)=FI 

IALAST=IA 

DO 1235 I=1917 

DO 1235 Jz1,IAMAX 

DO 1235 ΚΞ]1 917 

GAMMALI»J»K)=0.0 

THIS PT THE LANDING TRANSITION PROBABILITIES ARE COMPUTED. 
DO 300 IAA=1»1IAMAX | | 
DO 301 MI=1,IAMAX 

ΙΕΓ(ΙΑΑ-ΜΙ} 319432 933 

PRFMA(IAA»MI)=O. 


moo TO 301 х 


ММ1-МІ-1 
PRFMA(CIAASMI)ZP3MMI 
GO TO 301 

IAM1=TAA-1 


-MM1=MI-1 


BC(1)=1.0 


—PROD=FLOATF(IAA-MI). 


DO 50 IPz2,MI 
AIP=FLOATF(IP-1) 

PROD = PROD 41.0 

BC( IP) =PROD*BC(IP-1)/AIP 
IGO=IAA-MI 
PRFMA(TAA,»MI)= 
CONTINUE 
CONTINUE | 
THIS PT THE MAINTENANCE TRANSITION PROBABILITIES ARE 
DO 100 IT2152 


(BC(MI)*(1«e0-P)s**( IGO) ) *PX*MM] 


COMPUTED. 


43 


— 


0042 
0043 
0044 
0045 
0046 
0047 
0048 
0049 
0050 
0051 
0052 
0053 
0054 
0055 
0056 
0057 
0058 
0059 
0060 
0061 
0062 
0063 
0064 
0065 
0066 
0067 
0068 
0069 
0070 
0071 
0072 
0073 
0074. 
00751 
0076 
0077 
0078 
0079 
0080 
0081 
0082 
0083 


, 0084 


0085- 
0086 
0087 
0088 
0089 
0090 
0091 
0092 
0093 
0695 
0095 
0096 
0097 


25 
26 
28 
129 
1999 
14 


19 


10 


24 


20 
201 
102 
ШО 
100 


11 
T21 
17 


120 
AT 


86 


Γι τσ ο 25028 

TAU = TIME f 

GO TO 28 

TAU = UNITT-TIME 

υοιοι 131,0 

ОО 102 IJ=1,N 

ТЕЕ ОА 00.17 

ТЕ ИКО) 19,19,1999 | 
PTRM( Is TJ9IT) =EXPF(-FLAM*TAUXD) 
“στο 102 

PTRM(IsIJUsIT)=0. 

CONTO 102 

FIMISFLOATF(1J=1) 
PTRM(CISIJSIT)SEXPF (FLAM*TAU*F JM] ) 


GO TO 102 
О-о 
БО (ү =1„0 


PROD-FLOATF(I-IJ) 
DO 10 IP 521] 

AIP=FLOATF(.IP=-1) 

PROD = PROD + 1,0 

BC(IP) =PROD*BC(IP-1)/AIP 
ELT=EXPF(-FLAM*TAU) 

PTRM CIjsIJs IT)ZBCCIJ) (1, -ELT )*( I -IJ) *ELT (IJ-] ) 
GO TO 102 

IF(IJ-1-ID) 22,24,24 

CONTINUE 

CALL PID(ISIJSs IT) 

СОТО 102 


D=FLOATF(ID) 


ELDT=EXPF (-D*FLAM*TAU) 


FACT - 1.0 
A(1)21,0 
ΜΜ 5Ξ I-IJ 


DO 20 M=2,MM 
FACT=FACT+1.O 
A(M)=A(M=1)*FACT 
PTRMCISIJSIT)Z(D*FLAMXTAUJ*(I-IJ)*ELDT/ACI-IJ) 
CONTINUE 
CONTINUE 
CONTINUE 
GO τοι 
CONTINUE 
IF(IA-IALAST) О, 1215120 
IF (NEWN-NLAST) 11151175111 
IALPHZIALASTM1 $ GO TO 801 
CONTINUE 
THIS POINT THE LAUNCHING TRANSITION: PROBABILITIES ARE COMPUTED 
ОО 204 ІҒЕ-1»М 
IGM = XMINOF (IAsIFF) 
DO 203 IGz1,1GM 
IGM1sIG-1 
DO 202 IHs1jN 
IHM1=IH-1 
BPROD= ( (1. —PGAM)**]GM] )%*(PGAM**]HM]) 
IF(IG-IA) 91587584 
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“0131 


013 
013] 
013: 
013: 
01% 
0137 
0136 
013‹ 
0140. 
014; 
014% 
014: 
0144 
0148 
0146 
0141 
0146 
0145 
015€ 
0161 
0152 
0153 


— - е —— ——MÀ——  - - — - 00 





91 IFCIG+IHM1-IFF) 84582584 0154 
87 IFC IG+IHM1-IFF) 85585584 0155 
84 САММА(ІЕЕ» ІС» ІН) =0, 0156 

GO TO 202 0157 

82 BC(1)21.40 0158 
" PRODSFLOATF( IFF-IG) 0159 

DO 30 IPz2,IG 0160 
AIP=FLOATF(IP-1) 0161 
PROD = PROD + 1.0 0162 

-30 BC(TP)=PROD * BC(IP-1)/AIP 0163 
IHM1=IH-1 0164 
TEMPs PGAM**IHMI 0165 
TEMP1s(1«4-PGAM)**IGM] 0166 
BPROD = TEMP*TEMP1 0167 
САММА ( ІРЕ» 16» Ін) =ВС( ІС) *ВРВОО 0168 

GO TO 202 0169 

85 ЕВС(1)-1.0 0170 
PROD=FLOATF(IGMI-1) 0171 

DO 40 IP=2,IH | 0172 

AIP = FLOATF(IP-1) 0173 
PROD = PROD +1.0 0174 

40 FBC(IP)sPROD*FBC(IP-1)/AIP 0175 
САММА ( ІЕЕ» 16» ІН) =ЕВС( ІН) *ВРВОР 0176 

202 CONTINUE 0177 
203 CONTINUE 0178 
204 CONTINUE 0179 
REMOVE CARDS FROM HERE TO NO 999 IF PRINT OUT NOT DESIRED 0180 
PRINT 9>(((I>IJ>ITIPTRM(I>IJ>IT)>IT=192)>IJ=19N)>I=19NM) 0181 

9 FORMAT (1H1/(2(6H PTRM(I2>1H>12>1H>12>»3H) = El4+5))) 0182 
PRINT 99s t(LIFFsIGs IH»sGAMMALIFF»IG» IH) sIFF=1»N)s IG=1s IA) s IH=1»N) 0183 

99 FORMAT(1H1/(2(7H GAMMA(1I2>1H>I12>»1H>12>»3H) = El4+5))) 0184 
PRINT 999,(( IAASMISPRFMACIAASMIOS IAAZ15 IAMAX) $MIZ15 IAMAX) 0185 

999 FORMAT(1H1/(2(7H PRFMA(CI251H5I253H) = El4+5))) 0186 
NOW THE TRANSITION MATRIX MUST BF ZEROED 0187 
111 CONTINUE | 0188 
ОО 899 Ј=1,119 0189 

ОО 899 К=1›7 0190 

| 899 L=1517 0191 

899 PR(JsKs»L)=0.0 0192 
START COMPUTING THE ELEMENTS OF EACH ROWs I=IX+ (ALPHA - 1) X TOTAL- A/C 0193 

E DO 1000 IALPH=1> IALAST | 0194 

801 CONTINUE 0195 

DO 1100 IX=1»NLAST 0196 
COMPUTE THE P ELEMENTS OF THE IAPHsIX-ROW AND SUM THE ROW | 0197 
TSUM=0D0 0198 
I=IX+(TALPH-1)%N | 0199 
DO 800 IBETA=1>»IA : . 0200 
RSUM=0.0 0201 
DO 900 IY=1,NEWN 0202 
PR(I>IBETA>+IY)=008 0203 
ILIM=NEWN-IALPH-IX+2 | 0204 
PSUM=0.0 | 0205 
SUM#0e0 0206 
SUML=0+0 | 0207 
DO 500 IL=1sILIM | 0208 


KLIM=IX+IL-1 | 0209 
45 | 


701 
702 


700 


600 
500 


900 
800 


888 
L100 
1000 


889 
890 


898 
C NOW 
805 
807 


FORMAT ( 
CONTINUE 


See. | 
SUMM=0, 

DO 600 MI=1,IALPH 

SUMK=0. 

DO 700 IK=1s>KLIM 

IKPMI = IK +MI-1 

IF(IXPIL-NEWN) 70157015700 

IF(IKPMI-NEWN) 702,7025700 

GAMH=GAMMA ( ILIM> IBETA>IL) 

PTRMH1 = PTRM(IXPIL>IK»1) 

PREMAH = PREMA(IALPH»MI) 

PTRMH2 = PTRM(IKPMI>,IY>2) 

SUM = GAMH * PTRMHl * PRFMAH * PTRMH2 

SUMK =SUMK+SUM 

PSUM=PSUM+SUM 

CONTINUE 

SUMM = SUMM + SUMK 

CONTINUE 

SUML = SUML +SUMM 

PSUM2 =SUML р 
CONTINUE 

RSUM=R SUM+P SUM 

PR{IsIBETA» IY)=PSUM 

CONTINUE 
TSUM=TSUM+RSUM 
CONTINUE 
PRINT 888 » 


H Hf | 


TSUM y IALPHs IX 
RES UM = E 15.832 15.) 
CONTINUE 
REMOVE CARD FROM HERE TO 889 IF P MATRIX PRINT OUT NOT DESIRED 
DO 889 Jz1,17 
DO 889 К=1,7 
00889 1=1,17 
I=J+(K-1)*#N 
PRINT 8690, C PRSUTSLPSL)SLPZ1,.IAMAX) Ks Js L 
FORMAT (7E14.592HJ=12>5HK=1sA>2HL=12) 
ОО 898 1=1,119 
QJ(1I)=0.0 
MULTIPLY QI AND P TO GET GJ 
PRINT 807sKT»IALAST>»IA 
FORMAT(1H1»13HQ@ VECTOR CASE I3/// [5,15) 
DO 802 IBETA=1,7 
DO 902 IY=1,17 


SAT THIS POINT CALCULATE THE (IBETA,IY)TH ELEMENT OF THE QJ VECTOR 


2201 
2001 


J=IY+(IBETA-1)%N 
QP1=0. 

ОР=0. 

DO 2001 IALPH=197 
DO 2201 ТХИ 
I=IX+( TALPH-1)#N 
QP1=QI(I)#PR(I,IBETA»IY) 
QP=QP+OP1 

CONTINUE 

CONTINUE 

QJ(J)zQP 
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021 
021% 
021% 
021% 
021% 
021: 
0211: 
0214 
02 1% ή 
0221 
0228 


0221 
0221 
022: 
0226 
0221 
022€ 
022% 
023( 
0231 
0232 
02331 
0234 
02351 
0236 
0237! 
0238) 
0239 
02401 
02411 
0242] 
0243] 
0244! 
0245 
0246 
024T 
02481 
0249 
02501 
02511 
0252 
0253 
0254 
0255 
0256 
0257 
0258 
0259 
0260 
0261 
0262 
0263 
0264. 
02651 


PRINT 8882sIBETAsIYsJsQP 0266 





882 FORMAT(2I4s4H QJ(1393H )= El4,8) 0267 
B90? CONTINUE 0268 
$802 CONTINUE 0269 
E CHECK THE SUM OF THE Q VECTOR “ 70270 

QSUM=O. 0271 

DO 808 51.119 0272 

808 QSUM=QJ(U)+QSUM | 0273 
PRINT 8883,Q0SUM 0274 

883 FORMAT(6H QSUM= Ε15.9) 0275 
DO 333 I = 189119 | 0276 

ШЕ (І1-1)/17 0277 
FK=FLOATF(K) 0278 

x FMEAN= FK*QJ(I)+FMEAN | 0279 
333 CONTINUE | 0280 
TFLCsFMEAN*FLOATF( JCYCLE) | 0281 
PRINT 335,FMEAN 0282 

335 FORMAT( 17HMEAN A/C FLYING = F10+4) 0283 
STATEMENTS FROM THIS POINT TO THE CALL DRAW STATEMENT REFER TO GRAPH 0284 
JT(2)=ICODE(FMEAN) 0285 
JT(4)=ICODE(D) 0286 
FKT=FLOATF(KT) | 0287 
JT(6)=ICODE(FKT) 0288 
FNsFLOATF (NEWN-1) 0289 
JTC10)2ICODE(FN) 0290 
FIAASFLOATF(IA-1) 0291 
JTt12)5ICODE(FIAA) | қ 0292 

CALL DRAW(119»FJPLOT9QJ909094H J]1T909090909090N9898»909LAST ) 0293 
ЕМЕАМ - 0. | 0294: 
NEXT WE MUST MULTIPLY QJ AND P TO GET QK AND SO ON, səs (QK+,. +N) 0295 
| KT=KT+1 0296 
IF(KT-ICYCLE) 80358035806 0297 

ВОЗ DO 804 I=1,119 0298 

moe QO1(1)=QJ(1)} 0299 
| IALASTMIZIALAST | 0300 
| IALAST=IA 0301 

NLAST=NEWN 0302 

x GO TO 809 0303 
᾿ 806 STOP 06 | 0304 
| END ` 0305 
| SUBROUTINE PID(ISJSIT) 0306 
| COMMON FLAM,TIME 0307 
x COMMON PTRM, GAMMA PR I PRFMA ID 0308 
TYPE DOUBLE RC,RDC,PROD »DID3,NID4,DID5,DEXP 0309 

TYPE DOUBLE DAN »DID1 sDID2sSUM,DNsANM1 9 FAC sCOFsPSUMsPTRsFLAMs TAU sD 0310 
DIMENSION PTRM(1751752)5BC(11)5,5BDC(11) 0311 
DIMENSION GAMMA(1757917)s PREMA(797)s PR(11997517) 0312 
D=FLOATF(ID) | 0313 
IDP1=ID+1 0314 
IF(IT-1)25,25»26 0315 

25 TAU = TIME 0316 
GO TO 28 0317 

26 TAU= 1.-TIME 0318 
28 CONTINUE 0319 
IMDP1=I-ID 0320 
PTR=0.0 0321 
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=r = 


C 


CV CY CY CY £Y Су СУ СҮ СУ Су СУ СҮ СУ 


о. w ede $- - — ἀμ rl зше "πρ OA DI s... nn == — Ен. 
nn жаы - 


PSUM=0, 
DO 200 Νου τπτ ЖО 


DEVELOP N TAKEN J AT A TIME AND D TAKEN N AT A TIME 


10 


20 


a) 


207 


200 
102 


102 
101 


BC(1)=1.0 
PROD-FLOATF(NJ-J) 

DO 10 IPz2,J 
AIP=FLOATF(IP-1) 

PROD = PROD+1.O 
BC(IP)=PROD* BC(IP-1)/AIP 
PROD=FLOATF (IDP1-NJ) ` 

DO 20 IQ=2sNJ 
AIQ=FLOATF(IQ-1) 
PROD=PROD+1.0 

BDC( IQ) =PROD*BDC(IQ-1)/AIQ 
COF=BC(J)*BDC(NJ)%*(—1)%*(NJ-—-J) 

ANM1=FLOATF(NJ-1) 
DAN=D/(D-ANM]) 

DID4=DEXP (-FLAMXTAUXANMI]) 
0101 = ( БАМХ* ( 1-10Р1)) *0104 
SUM=06 

DN=06 

DO 201 K=1>IMDP1 

FAC=l. 

ΚΜΙΞΚ-1 

PROD=0. 

DO 11 IK=1»KMl 
PROD=PROD+l. 

FAC=FAC*PROD 

IMIDK=I-ID=K 

SUM= ( (FLAM#D*TAU)##KM] )#DÄN**IMIDK / FAC 
DN=DN+SUM 

DID3=DEXP (-FLAMXD*TAU) 
DID2=0N*DID3 
DID5-DIDI-DID2 
PSUM=COF*DIDS 

PIR πο πο ο 

CONTINUE 

PTRM(1»J»IT)=PTR 

CONTINUE | 

CONTINUE 

END 

FUNCTION KRAN(As$B,. IY) 


THIS ROUTINE RETURNS AN UNIFORMLY DISTRIBUTED RANDOM INTEGER 


THIS ROUTINE RETURNS A INTEGER RANDOM NUMBER’ .GE. TO A 

“Т, В 

A = BOTTOM LIMIT (INCLUDED) FOR THE RANDOM NUMBER 

В = TOP LIMIT (NOT INCLUDED) FOR THE RANDOM NUMBER 

SET IY ONLY ONCE IN MAIN PROGRAM FOR EACH SET OF RANDOM NUMBERS 


SOME GOOD STARTING VALUES FOR IY FOLLOW 
13421773 


33554433 
8426219 
42758321 
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26237485 
62104023 | 
ANY OF THESE MAY BE USED 


ВЕТ ROUTINE MAY BE USED IN FORTRAN 60 OR 63 


ozo IY 

IY = IY - (1Y/67108864) * 67108864 
FY = IY 

- KRAN = FY/67108864. * (B-A) + A 


i FINIS 
XECUTER. 


0378 
0379 
0380 
0381 
0382 
0383 
0384 
0385 
0386 
0387 
0388 
0389 
0390 
0391 


— — - = ` ai atia PPP o == 


Q VECTOR CASE 1 


H 
| * 





5 5 | 
1 ТТЫ 20970217306 08 4 1 0/0 52 )-5,7976429ЭЕ-Ц 
1 20/0 2 )-3.020208351Е-02 4 2 0J( 53 )=3.03224061E-02 
1 à QJC 3 )=1.72436823E-02 4 3 0J( 54 )=1.16230913E-02 
1 4 0J( 4 )z6.68152510E-03 4 4 QJ( 55 )=2.87808733E-03 
i 5 QJ( 5 )=2.02525670E-03 4 5 0./( 56 )=5.26888829Е-04 
1 6 QUC 6 2-5.29756823Е-04 4 6 QJ( 57 )=7.84720454E-05 
1 7 0Ш8( 7 )z1.30002033E-04 4 7 0J( 58 )=1.04516559E-05 
al 8 QJC 8 ):=3.22530060E-05 4 8 QUC 59 )=1.39145026E-06 
1 9 0./( 9 )=8.58857000Е-06 4 9 0./( 60 )=2.08960624Е-07 
1 10 0J( 10 )=2.86285667E-06 4 10 QJ( 64 2:4.55354654Е-08 
1 11 QJ( 11 ):8.34999861E-07 4 11 QJ( 62 ):7.2555772 98 209 
1 12 0J( 12 )=2.08749965E-07 4 12 0.( 63 )=9.06947154Е-10 
1 13 QUC 13 )=4.34895761E-98 4 18 QJ( 64 )=7.55789295Е-11 
1 14 QJ( 14 )=7.24826263E-09 4 14 QJ( 65 )=3.14912206E-12 
1 15 QJ( 15 )=9.06032836E-10 4 15 0J( ος πε 0 
1 16 QJ( 16 )z7.55027363E-11 4 16 QJ( 627): 0 
1 17 QJ( 17 )z3.14594735E-12 4 47 QJ( 68 )= 0 
2 1 0/( 18 )=5.02519433E-02 5 1 QUC 69 )=4.9164283E-01 
2 2 QJC 19 )=3.07381552E-02 5 2 0./( 70 251.45718127Е-01 
2 š QJ( 29 )=1.54385405E-02 5 3 QJ( 71 )=2.69531853E-02 
2 4 0J( 21 )=5.18865107E-03 5 4 QJC 72 )=3.7879962E-03 
2 5 QJ( 22 )=1.34374437E-03 5 5 QJ( 73 )=4.23732653E-04 
2 6 QJ( 23 )=2.96164302E-04 5 6 QJC 74 )=3.99381550E-05 
2 7 QJ( 24 )=z6.07227021F-05 5 7 QJ( 75 )=3.41914062E-06 
2 8 QJ( 25 )=z1.26203649E-05 5 8 QJ( 76 )=2.98561135E-07 
2 9 QJ( 26 )=2.86289508E-06 5 9 QJt 77 )=3.11088104E-08 
2 10 0J( 27 )=8.35011051E-07 5 10 QUC 78 )=4.95114753E-09 
2 11 0J( 28 )=2.08752763E-07 5 11 QJ( 79 )=5.92108046E-10 
2 12 QUC 29 )=4.34901589E-08 5 12 0J( 80 )=4.72886874E-11 
2 13 QJ( 30 )=7.24835981E-09 5 13 0J( 81 )=1.89139810E-12 
2 14 QJt 31 )=9.06044977E-10 ET οὓς 82 - 0 
2 15 0J( 32 20:7.55037481E-11 5 15 QJ( 83 )= 0 
2 16 QJ( 33 )=3.14598950E-12 5 16 QJ( 84 )= 0 
ο. 54 )- 0 ο OJ ВБ ЕЕ 0 
3 1 0J( 35 )=3.40989043E-02 6 1 QJC 86 )= 0 
3 2 QJC 36 )=3.07597531E-02 6 2 QJ( 87 у= 0 
3 šÉ QJ( 37 ):1.35265107E-02 6 $ 0/С 89 oe 0 
3 4 0J( 38 )=3.91531781E-03 5 4 QJ( 89 )= 0 
É 5 0Ј( 39 )=8.57092801Е-04 6 9 οί 90 ἘῈ 0 
5 6 QJ( 40 )=1.56625199E-94 6 6 OIC WIE 0. 
3 7 QJC 44 )=2.62271034E-05 6 το 9275 0 
3 8 QJt 42 )=4.44195168E-96 6 8 0J( 93 De 0 
3 0290901545 )28.3512595 9E 707 6 9 QJC 94 )z 0 
3 10 0J( 44 )=2.08781489E-07 6 10 QJ( 95 >= 0 
3 11 QJt 45 )=4,349614353E-08 6 11 0J( 96 )= 0 
š 12 QJ( 46 )=7.24935723E-09 6 12 QJ( 97 )= 0 
3 13 QJ( 47 )=9.906169656E-19 6 13 0J( 98 )=- 0 
3 14 QJ( 4B )Ξ7.55141380Ε-11 ' 6 14 QJ( 99 )= 0 
š 15 0J( 49 )=3.14642242E-12 6 15 QJ(100 )š 0 
9 ТӨШ "ορ. 0 6. 16 ο που }- 0 
ο το Eg E 0 A RA 0 
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103 3 
ο. (104 ) 
οὐ (105 )= 
ου (406 = 
QJ(107 )= 
0J(108 ) 
0J(109 ) 
000110 

ο (1. 

τ ΠΠ 6 

MOJ ILS 

no GJ6]479)- 

dS up GE Ба дан 

14 QJ(116 )= 

1> 07M 7 )- 

16 0J(118 ) = 

17 QJ(119 )- 
OSUM=1.000000000E 00 

EAN A/C FLYING 53.1666Е 00 


OO мо хл οἱ Νε"... 


Y Y Y Y Y Y Y Y Y Y Y Y Y Y Y YY о ы 
OO O OO OO OO о о о о осо 


GRAPH TITLED | 
ECA/C) =3.17E+00SPOTS 
J VS QJ VECIOR N 
HAS BEEN PLOTTED. 


1.00E*00 
4.00E*00 


.00Е%00 Т 
.60Е%01 А 


и 0 
r © 
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АРЕЕМОТХ ГУ 


SAMPLE RESULTS 


The following pages present the values of the elements of the proba- 
bility distribution vector (QJ) and its graphical plot for five consecutive 
iterations, i.e., Qx EN for n = 1, 2,,,.,5. The inputs are those 
shown on the first page of Appendix III between statement No. 30 and 
No. 31. The printouts of the transition matrices and their computational 
elements are omitted. The plot was made using the DRAW subroutine in 
the U.S. Naval Postgraduate School computer facility library. Each 
vector printout contains the values of all 119 states possible (7 x 17) and 
is headed by the past value of A+ 1 and the next value of A+ 1. The two 
indices preceding each element represent B + l andj + l, in the notation 
of section 3. For example, in the first row on the next page, the "] 1" 
indicates that the probability of being in state (0, 0) after one iteration 
is 2 ,026, where the value of A is 4 over the first iteration. Each graph 
is labeled with the expected value of a/c flying, the number of mainte- 
nance spots available, the vector number (T), total number of a/c avail- 
able (N), and the desired number of a/c on station (A). The "E" notation 
indicates the power of 10 to multiply by. This sample run demonstrates 


the loss in total a/c and variable a/c on station. 
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_ 995 _ 


834 005 


803. 


002 


. 
—O  — o ——— s — nn nn nn m - — —  — —— —] — - 


TX ПЛ x 


0090 002 984 096 9096 010 019 


909 


M"SCRLE - 2.00E4-91 UNITS^INCH. 
Y-SCALE = 1.002 -91 UNITS-INCH, 


Е(В/С) -3.12Е5--005Р015 


= &.00r 4-00 Т 
J YS QU VECTOR E. 


1.00E+0 1 


"Ее. 


D y 

е 

E- O 
ОН о 


anteca mapaa 


Q VECTOR CASE 2 





5 5 
1 19005 1 )=9.12184624Е-06 4 4 00( 52 ):7.25055242Е ПО 
1 2 QJ( 2 )=7,94491028E-06 4 2 QJ( 53 )524.45693258E-03 
1 3 0J( 3 >)=3.41380533E-06 4 à QJ( 54 )=1.28968676Е-03 
1 4 QJ( 4 )=9.83913729E-07 4 4 0J( 55 )=2.35500532E-04 
1 5 0./( 5 2:2.20530359Е-07 4 5 QJ( 56 )=5.11156232Е-05 
1 -6 QJ( 6 )=4.28567319E-08 4 6 QJC 57 )=3.30021596E-06 
1 7 090 7 3:7.97994342Е-09 4 7 QUC 58 155.15123146Е-ТПИ 
j 8 QJC 8 )=1.56348438E-09 4 8 QJ( 59 )=3.13674353E-08 
Е 9 QU( 9 )=3,46529220E-10 4 9 QJC 60 )=3.81044171E-09 
1 10 0J( 10 )=1.01221242Е-10 4 10 0J( 61 )=6.98663591E-10 
101 0001100 2 δες SUSE 0] 4 11 QJ( 62 )=1.02568738E-10 
1 12 OJ 12 )=5.63385201Е-12 4 12 0306.63 У=1.13092255Е=1 
fe 15 QUC 1343 )21 122986 БЕ 222 4 418 QJ( 64 ):8.32962386E-13 
1 14 QJ 14 )=1.48426307E-13 4 14 QJC( 65 323.9 755 9,5 9 SE 
10315 DJ 1 1.0141 9 BE me Ane 4 5 ος. EE 0 
1 16 QJ( 16 )z1.17030587E-15 4 16 QJ( 67 = 0 
1 17 QJ( 17 ):=4.24536876E-17 4 17 QJ( 68 )= 0 
2 1 0J( 18 )=1.69105396Е-04 5 1 0./( 69 2:7.70076632Е-01 
2 2 QJ( 19 )=1.31958642E-94 5 2 QJ( 70 )=1.88159529E-01 
2 3 QJ( 2g )=4.93425071E-05 5 š QJ( 71 )=2.32480412E-02 
2 4 QJ( 21 )=1.22416819E-05 5 4 QJ( 72 )=1.93788734E-03 
2 ο ο, ορ το 3169525097 бб "5 5 0./( 73 )=1.22762068Е-04 
2 6 QJ( 28 )=3.73462197E£-07 5 6 0Ј( 74 )=6.39420794Е-06 . 
2 7 Qul 24 )-5.70661929Е-08 5 7090 75 )=3,02738808E-07 | 
2 8 QJ( 25 )=9.22154646E-99 5 8 0./( 76 )=1.56312395Е-08 
2 9 QJ( 26 20-1.72702939E-09 5 9 QJ( 77 )=1.13519205E-09 
2 10 0J( 27 )=4.41481881E-10 5 10 0J( 78 )=1.50627061E-10 
2 11 010 28 )59,65869676Е-11 5 11 0J( 79 )=1.51200491E-11 
2 12 0J( 29 )=1.75883018E-11 5 12 0J( 80 )=1.02093434E-12 
2 13 QJ( 30 )=2.56023495E-12 5 13 0J( 81 ):3.47880865E-14 
2 14 QJ( $1 )=2.79434001E-13 5 14 QJ 82E 0 
24 15 QJ( 32 )=2.033930353E-14 5 15 QJ( 83 )= 0 
2 16 QJ( 33 )=7.41007441E-16 Ὁ ο ος 84 у= 0 
елек Say 0 DE "170908559" 0 
3 1 0J( 35 )=1.42525319E-03 6 1 0J( 86 )= 0 
3 2 0J( 36 )=9.83757697Е-04 6 ους 5 15 0 
3 5 0Ј( 37 )=3.24331350€-04 6 ο ος 88 E 0 
3 4 QJ( 38 )=6.90234010E-05 ó 4 QJ( 89 )= 0 
3 5 QJ( 39 )=1.09403657E-05 6 DESC 99 у= 0 
3 6 0.( 40 )=1.44032446Е-06 6 6 QJ( 91 )= 0 
É 7 QJ( 41 )=1.76353392E-07 6 EU 92 E 0 
3 8 0 442 )-2.28151763Е-08 6 8 οὐ 93 )= 0 
3 9 QJ( 43 )=3.51078290E-09 6 9 0J( 94 Ye 0 
3 10 QJ( 44 )z7.70559333E-10 ern 0) 95 D 0 
š 11 QJ( 45 )z1.40851171E-10 6 11 02 ο. - 0 
3 12 QJ( 46 ):2.05928494Е-11 6912 0J 9795€ 0 
3 13 GUC 47 )=2.25869291E-12 © 6 13 0J( 985% 0 
3 14 QJ( 48 3:1.65302266Е-15 6 14 QJ( 99 )= 0 
3 15 QJ( 49 )z6.05799194E-15 6 15 QJ(100 )š 0 
5-16 0/5 51 )= 0 6 16 0J(101 )= 0 
uy JC )- 0 6 17 QJ(102 )- 0 
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ою о JO MU + GÇ (fÑ > 


οὐ (103 
090104 
ο (105 
ο (106 
091107 
QJ(108 
QJC109 
lin 
11 
QI TI? 
061153 
QJ(114 
MOTO 
QJ(116 
UG 10387 
QJ(118 
QJC119 
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) = 


)* 
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у 


jum 
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OSUM=1.000000000E 00 


EAN A/C FLYING =3.9799E 00 
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225 
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892 


001 
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000 002 004 005 
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ὃ. 


QQE 4-00 
1.66Е+01 





A e - sve tma . + —— = ae . — δα AX PP m s ww — a — —sa c s ar 
ыы ғ "ч. IR FE Fe En c ym miam pr. o y "алад ттер” «еф тұры ар ен tec ^ 
5 . ` e < м ад Ў A eons 


Q VECTOR CASE 3 


5 7 
1 ШЕЛІ 1 2-2 8556821 46-15 4 ο 5 )=1.66499369E-02 
1 2 QJ( 2 ):2.15851583Е-09 4 ο στο ):9.01025249Е-035 
1 0 3 )=7.94343376Е-06 4 3 QJ( 54 )=2,23842140E-03 
1 4 QJ( 4 )=1.91420746Е-05 4 4 QJ( 55 )=3.40006575E-04 
ik 2 QJ( 5 )=3.49125464E-07 - 5 QJ( 56 )=3.57140995E-03 
i 6 QJ( 6 )=9,37405166Е-08 4 6 QJ( 57 )=2.82593562E-05 
1 7 020 7 3:7.77939446Е-09 + 7 QJ( 58 )=1.85894467Е-07 
1 8 00( 8 3:1.187566095Е-02 4 8 QJ( 59 )=1.193555089E-08 
1 9 QJ( 9 )=2.11283841E-10 , 5 9 QJ( 60 2029.515295527E-10 
10 0.0 10 )229.18992660E-11 * 4 10 0J( 61 )=1.28538145E-10 
1 11 0J( 11 )=1.08528620E-11 4 11 QJ( 62 )=1.28940845E -11 
Ie DE 12 )=1.8773930%E-1?2 4 12 QJ( 63 )=8.53320924E-13 
1 135 QJ( 13 ):=2.57762883E-13 4 13 0J( 64 )=2.79262670E-14 
1 14 0J( 14 )=2.63176291E-14 4 14 QJ¢ 65 35 0 
1 15 0J( 15 )=1.77499822E-15 4 15 QJ( 66 )= 0 
1 16 QJ( 16 ):5.92695364Е-17 4 16 QJ( 67 )= 0 
EJ 7030 279 0 4 17 QJ( 68 )=: 0 
2 1 0J( 18 )=4.81711691E-04 2 1 QJt 69 )24.98253597E-02 
2 2 0J( 19 )=9.26173533E-94 2 2 QJ( 70 )=2.42738832E-02 
2 $ 94( 20 )=1.04905695Е-04 5 3 0./( 71 ):5.55895703Е-03 
2 4 QJt 21 )=2.16142339E-05 5 4 QJ( 72 2:7.10282975Е-04 
2 5 QJ( 22 )=3.28306882E-0% 5 5 0J( 73 )=6.33588530E-05 
2 6 QJ( 23 )=4.09033643E-07 5 6 0./( 74 )=4.082328683Е-06 
2 7 0)( 24 )-4.68083562Е-03 5. 7 0J( 75 )=2.05024689£-07 
2 8 QJ( 25 )=52,62586375Е-09 5 8 QJ( 76 2:9.26759069Е-09 
2 9 QJ( 26 )-8.1ῃ258169Ε-1ῃ 5 901( 77 )=5.04461339Е-10 
2 10 0J( 27 )=1.69317736E-10 5 10 0J( 78 )=5.13281492E-11 
2 11 QJ( 28 )=2.92764115E-11 5 11 QJ( 79 )=3,44456222E-12 
2 12 0J( 29 )=4.01871703E-12 5 12 QJ( 80 )21.14275058E-15 
2 13 0J( 30 )=4.10297009E-13 5 13 QJ( 81 )=: 0 
2.14 QJ( $4 )=2.76746676E-14 5 14 QJ( 82 )= 0 
2 15 QJ( 32 )=9.24174352E-16 5 ους = 0 
PNE 09006550) = 0 5 16 QJ( 84 ) = 0 
а са Ü 5> 17 QJ( 85 E 0 
5 1 QJ( 35 )=3.67589953E-03 6 1 QJ( 86 )=1.0430002E-01 
3 2 QJ( 36 )=2.22029313E-03 6 2 QJ( 87 )=4.52895010E-02 
3 3 QJC 37 )=6.25702594E-04 6 3 QJ( 88 )-8.79453771Е-03 
3 4 0J( 38 )=1.10292583E-04 6 4 090.( 89 )=1.006905238Е-03 
3 5 QJC 39 2:1.38900387E-05 6 5 0Ј( 90 )=7.55512291Е-05 
3 | 6 QJ( 49 )=1.38024736E-06 6 6 QJ( 91 )=3.92418576E-06 
3 7 QJ( 41 )=1.21057358E-07 6 7 QJ( 92 )=1.47699319E-07 
3 8 0J( 42 )=1.09262020E-08 6 8 QJ( 93 )=4.43630448E-09 
3 9 QJ( 43 )=1.21362937E-09 6 9 QJ( 94 )=21.43899935E-10 
$ 10 QJ( 44 )=2.09057619E-10 6 10 Οὐ 95 ):9.80396467Е-12 
ὁ 11 0J( 45 )=2.860955733E-11 6 11 QJ( 96 )=3.3(144042E-13 
12 QJ( 46 )=2.91366712E-12 6 12 QJ( 97 )= 0 
ὁ 13 QJ( 47 )=1.96282619E-13 6 13 0J( 98 )= 0 
ὁ 14 QJ( 48 )z6.55549374E-15 6 14 0J( 99 )= 0 
S 15 QJ( 49 )= Ü 6 15 QJ(100 ) = 0 
БЕС OJ( 59 )= Ü 6 16 QJ(101 )= 0 
EE QUIC 54 2- 0 6 17 0J(102 >= 0 
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7 1 GJ(108 ):5.55639840E-01 
7 2 QJC104 )=1.48858751E-01 
7 3 0J3(105 )=1.84227232E-02 
7 4 09(1106 )51.57794578Е-03 
7 5 QJC107 )=6.83380035E-05 
7 6 QJ(108 )=2.32262953E-06 
7 7 QJ(109 )=5.43332378E-08 
7 8 QJ(110 )=8.80653967E-10 
7 9 QJ(111 )=1.18536160E-11 
7 10 0J(112 )=3.82356527E-13 
7 een 0 
О ОЈ .π- 0 
715 ο αρ. - 0 
7 “14 IM NE - 0 
7 —15980J0] - 0 
7 16 ade 0 
 Ἱπται Е 0. 
OSUM=1.000000000E 00 
"EAN A/C FLYING =5.5636E 00 
GRAPH TITLED 
Е(А/С) =5.56E+Q0SPOTS = 8.00E+00 T = 3.00Е%00 
J VS QJ VECTOR N = 1.50E+01 AE MO EDO 


HAS BEEN PLOTTED. 


58 








GBS 


99+ 


023 


892 


221 


020 


009 


002 094 


A OCALE = 2.00E+01 UNITS INCH, 
Y-SCRLE = 1.006 -@1 UNITS~INCH. 


006 


E(Q^C) -ο.50ξ--ΘΘ5ΡΩ 19 
J US QJ VECTOR 


59 


008 


М 


8.90. +08 
= 1,08. +01 


019 


912 


| ——À --- 


A MÀ ———  — a 


Q VECTOR CASE 4 





7 57 
1 т 9-3. 5700636 0Е- ПМ 4 1 QJ( 52 )=8.03893389E-02 
1 2 QJ( 2 )=2.66337542Е-04 ` 4 2 QJ( 53 )=3.58842827E-02 
1 $ 0/0 3 )=7.61707863Е-0Э 4 3 0Ј( 54 )=7.24458723Е-03 
1 4 QJC 4 )=1.37453763Е-(05 4. 4 QJ( 55 )=8.75756513F-94 | 
1 5 QU 5 )=1.79617995Е-06 4 5 QJ( 56 )=7.99643835E-95 | 
1 6 040 6 3:1.89112432Е-07 4 6 0./( 57 )=4.12757543Е-06 | 
1 7 QJ( 7 )=1.81098891Е-08 4 7 QU 58 )=1.85478755Е-07 | 
1 8 QJ( 8 )=1.84326727E-09 4 8 QJ( 59 )=7.45697706E-09 i 
1 9 QJ( 9 )=2.33815659E-10 4 9 QJ( 69 )=23.67973791E-10 i 
4 109 9/0 10 72:4.51125756Е-11 73 "4 10 0/0 61 7-%956621692Е ul | 
1 11 Q0J( 11 )=7.18762832E-12 4 11 0J( 62 )=2,27683580E-12 
πιο ο 125-997 64267 Geos 4 12 QJ( 63 )=7.17866372E-14 BE 
1 15 @Ј( 13 )=8.51451858Е-14 2015 QJC 6355 0 
1 14 QJt 14 )=5.27326607E-15 4 14 070 65 95 0 | 
1 15 QJ( 15 )=1.61678538E-16 4 15 QJ( 66 )= 0 mE 
1 16 QJ 16 )= 0 4 2 80 (02 26 0 
1 817 DJ a hs 0 τι; OU 68 у- 0 
2 1 0J( 18 )=5.18126447E-04 5 1 0J( 69 )=6.25520845E-01 | 
2 2 GUC 19 )=2.86433504E-03 5 2 QJ( 70 )=1.72496943E-01 | 
2 3 QJ( 20 )=7.33497609E-04 5 3 QJ( 71 )=2.2094320JE-02 | 
2 4 QJ( 21 2-1.16540651E-04 5 4 QJ( 72 »)-:1.7$10455V E SIDES | 
2 5 Ш  22)51.40125851Е-05 5 50. 73 )=9.18785067Е-05 № 
2 б 0Ј( 235 )=1.15029797/Е-06 5 6 QJ( 74 )=5.46552930Е-06 
2 7 QJ( 24 228.54466567E-98 5 7 QJ( 75 )=9.61462191E-03 
2 8 0./( 25 )=6.65517506Е-09 5 8 QJ( 76 )=2.135377/50585309 
2 9 QJ( 26 )=6.57969504Е-10 5 9 01( 77 ›=5. $1168 753Е αι 
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